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We report on generic relations between fractional flow and pressure in steady two-phase flow in 
porous media. The main result is a differential equation for fractional flow as a function of phase 
saturation. We infer this result from two underlying observations of steady flow simulations in 
two and three dimensions using biperiodic boundary conditions. The resulting equation is solved 
generally, and the result is tested against simulations and experimental relative permeability results 
found in the literature. 



Two-phase flow in porous media is possible to classify 
as either steady or unsteady flow. While steady flow is 
the one that is statistically invariant in time, the typ- 
ical unsteady flow is the displacement of one phase by 
another. In particular, many authors have studied the 
motion of the displacement front. This was done by 
experiment network modeling Q], lattice Boltzmann 
methods[|| and also statistical models: diffusion limited 
aggregation [|] and invasion percolation j|. The study of 
steady flow is often associated with finding the relative 
permeability curves, although this concept can be de- 
fined and is used both for steady and unsteady flowJB], ||. 
Comparisons between relative permeability obtained by 
the two methods exist in the literature^. However, in 
the complexity of two-phase flow in porous media, pure 
displacements and true steady flow can be considered to 
be two limiting situations. Thus, by nature the two sit- 
uations are very different, and one should expect differ- 
ences in the relative permeability curves. Nevertheless, 
unsteady methods are generally preferred because steady 
methods take long time and are more expensive. The re- 
sults that we present here concern steady-state flow. A 
lattice-Boltzmann approach to steady-state flow is found 
inS. Further, a very interesting study using a network 
approach is found in|^]. 

The body of knowledge in the field is immense B [To] . 
Pore-network modeling is popular and a good overview 
of the state-of-art is given inO]. Basic mechanisms of 
displacement on pore level are largely known, and much 
qualitative knowledge exists about properties on larger 
scale. However, there is a good deal of work left to do in 
order to bridge the gap between these two worlds. Nu- 
merical simulations is a unique tool to that respect being 
very well controllable and more precise. In particular, 
we find that while experimental data points often are 
few and scattered, we have the possibility to generate 
sufficient data points to obtain good statistics. Further, 
in simulation series we can control the ensemble com- 
pletely. Commonly used ensembles for relative perme- 
ability or fractional flow curves as functions of satura- 
tion, are constant global pressure, constant total flux or 
constant capillary number(Ca). We argue that the choice 
of ensemble can be of major importance and is often un- 



derestimated. Our results are presented in the formalism 
of fractional flow and globally applied pressure. We find 
that when using the constant Ca ensemble, there exist 
a generic differential equation for the fractional flow. In 
turn a general solution is found. 

The results of the paper is based on a network simula- 
tor for immiscible two-phase flow. This line of modeling 
which is based on Washburn's equation |l2| dates back to 
the work of several groups in the mid 1980's|l| [l|, [§. 
Our model is a continuation of the model developed by 
Aker et aZ. jlfj, |l7|]. Although a thorough presentation 
can be found in Jl8| , we provide a brief resume of main 
aspects for clarity. 

The porous media are represented by networks of 
tubes, forming square lattices in two dimensions (2D) 
and cubic lattices in 3D, both tilted with respect to the 
imposed pressure gradient and thus the overall direction 
of flow. We refer to the the lattice points were four(2D) 
and six(3D) tubes meet as nodes. Volume in the model 
is contained in the tubes and not in the nodes, although 
effective node volumes are used in the modeling of trans- 
port through the nodes. For further details, see @. Ran- 
domness is incorporated by distorting the nodes on ran- 
dom within a circle(2D) and a sphcrc(3D) around their 
respective lattice positions. This gives a distribution of 
tube lengths in the system. Further, the radii are drawn 
from a flat distribution so that the radius of a given tube 
is r G (0.1Z, 0.4Z) where I is the length of that tube. 

The model is filled with two phases that flow within 
the system of tubes. The flow in each tube obeys the 
Washburn fl2f equation, q = — (ak/fi)(Ap — J^Pc)/^- 
With respect to momentum transfer these tubes are 
cylindrical with cross-sectional area a and length I. The 
permeability is k = r 2 /8 which is known for Hagen- 
Poiseulle flow. Further, ijl is the viscosity of the phase 
present in the tube. If both phases are present the vol- 
ume average of their viscosities is used. The volumet- 
ric flow rate is denoted by q and the pressure differ- 
ence between the ends of the tube by Ap. The sum- 
mation is the sum over all capillary pressures p c within 
the tube. With respect to capillary pressure the tubes 
are hour-glass shaped meaning that a meniscus at po- 
sition x £ (0, 1) in the tube has capillary pressure: 
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p c = (2-f/r)[l — cos (2irx/l)], where 7 is the interfacial 
tension between the two phases. This is a modified ver- 
sion of the Young-Laplace law(|, |i~7| ]. 

Biperiodic(2D) and triperiodic(3D) boundary condi- 
tions are used so that the flow is by construction steady 
flow. That is to say, the systems are closed so both phases 
retain their initial volume fractions, i.e. their satura- 
tions. These boundary conditions are in 2D equivalent 
to saying that the flow is restrained to be on the sur- 
face of a torus. The flow is driven by a globally applied 
pressure gradient. Usually invasion processes are driven 
by setting up a pressure fall between two borders, inlet 
and outlet. Since, by construction, the outlet is directly 
joined with the inlet in our system, we give instead a so- 
called global pressure drop when passing this line or cut 
through the system]!^, Effectively we put restraints 
on the pressure gradient that is experienced throughout 
the network. Integration of the pressure gradient along 
an arbitrary closed path one lap around the system, mak- 
ing sure to pass the 'inlet-outlet'-cut once, should add up 
to the same global pressure fall. 

For a given value of the global pressure fall, and given 
the instantaneous distribution of the phases in the sys- 
tem, the pressure field is calculated numerically, and the 
flow distribution at that moment is known. Given the 
flow field the whole system is updated according to the 
Euler scheme. All interfaces, menisci, which are the es- 
sential entities keeping track of the distribution of the 
phases, are moved according to the fluid velocity in the 
tube where they are situated. After having updated the 
hole system, the new distribution of the phases might 
give new effective viscosities in the tubes and different 
capillary pressures across the menisci, leading to a re- 
calculation of the coefficients in the set of equations for 
the pressure. The tricky part is how to model the trans- 
port of the menisci across the nodes, and this is done 
by applying a set of rules that assures volume conserva- 
tion of both phases and that are acceptably close to the 
real world. We refer to the previous study for details on 
this@. 

All simulations being the basis for our conclusions are 
done in series where the saturation is varied as the inde- 
pendent variable and the capillary number as well as the 
viscosity ratio M — [i nw /[i w is held fixed. We define the 
capillary number as Ca = QtotA*eff/(^7) where E is the 
cross-sectional area of the entire network, Qtot is the to- 
tal flux through the network and ^j, c g = fJ, nv ,F DW + fi w F w . 
Here, F nw and F w are nonwetting and wetting fractional 
flow respectively. In invasion studies one often employs 
either the wetting or the nonwetting viscosity and refer 
to the wetting or nonwetting capillary number, respec- 
tively. However, in steady flow it is appropriate to take 
the volume average. The resulting dependent variables 
are the fractional flow of the two phases and the globally 
applied pressure. These two variables are in a one-to- 
one relationship with the two relative permeabilities, but 
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FIG. 1: The figure shows that the relation between pressure 
and the derivative of the fractional is of the form in Eq. (Q). 
The curves are for viscosity matching phases for four different 
capillary numbers in 2D and two capillary numbers in 3D. The 
nonwetting saturation 5 n w serves as parameter for the curves, 
and the intervals in which the data is plotted are as follows; 
(2D) Ca = 3.2 x 1(T 2 : S nv> G (0.14, 0.93), Ga = 1.0 x 10" 2 : 
G (0.14, 0.83),Ca = 3.2 x 10~ 3 : S nw G (0.14, 0.83), Ga = 
1.0 x 10~ 3 : G (0.20,0.73), (3D) Ca = 1.0 x 10~ 2 : S nw G 
(0.10, 0.90), Ca = 1.0 x 10" 3 : SW G (0.20,0.75). The two 
curves in 3D are lowered by one unit for clarity. 



we will stick to the former notation here, for clarity and 
convenience. 

From the simulations two observations are made from 
which we infer the following differential equation for the 
fractional flow as a function of saturation; 
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(1) 



The saturation is the nonwetting saturation S — S nw and 
the fractional flow is the nonwetting fractional flow F — 
-Fnw(Snw) here and in illustrations unless where explicitly 
stated otherwise. However, the relation is equivalent in 
form for the wetting counterpart F W (S W ), only the values 
of a and b change. 

The two observations on which Eq. (|l|) is based are of 
a phenomenological character. The physical problem is 
intractable from an analytical point of view. By means 
of simulations the evolution and properties are monitored 
in detail. The first observation is that the derivative of 
the fractional flow is related to global pressure by 



(2) 



where P denotes normalized pressure, i.e., the actual 
pressure for the fixed flux at a given saturation divided by 
the pressure for single phase flow at the same flux. This 
result has been reported on in great detail for the case of 
viscosity matching phases) 20 . Those results were solely 



based on simulations in two dimensions. However, the 
result is extendable to 3D, see Fig.[l], as well as viscosity 
contrasts. 
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FIG. 2: The figure shows fitted curves to P(Kiw) according to 
Eq. (^|). The results are from the same series of simulations as 
the curves in 2D in Fig. [j], only showing results for two lower 
capillary numbers. The shown region of F nw is the one where 
the fit is good. For these capillary number this region is large, 
for the two higher values of Ca shown if Fig. |lj the region of 
good fit is; Ca = 3.2 x 10" 2 : F nv E (0.15, 0.45), Ca = 1.0 x 
10 -2 : -Fnw S (0.10, 0.50). The data is from five realizations of 
the porous network giving slightly shifted curves as we can see. 
The scattering of the data becomes larger for smaller capillary 
numbers due to increased hysteresis and history effects. 



The second observation is that for a broad range of 
saturation values the pressure forms a quadratic function 
of the fractional flow. In a mathematical formulation this 
result becomes 



P(F) = a'F 2 + b'F + c! 



(3) 



The range of validity of this equation is illustrated by 
the quadratic fits in Figs. || and [| In Fig. | the vis- 
cosity ratio is unity, and we observe how Eq. (||) holds 
for each fixed value of Ca. Likewise, we show in Fig. 
H quadratic fits to P(F)-curves for one given capillary 
number Ca = 3.2 x 10 -3 , varying the viscosity ratio from 
M = 30 to M = 1/30. The point is that for most of the 
range that we are interested in, it is sufficient to expand 
P(F) to second order in F around the maximum point. 
Alternatively, an expansion of P(S) in S turns out to re- 
quire higher order terms, and is hence less useful. Other 
combinations of capillary number and viscosity ratio have 
been performed and they are found to agree with this re- 
sult. The shifting of the curves and detailed dependence 
on Ca and M are to be discussed elsewhere pl|. 

Recall that S is the independent variable upon which 
both F and P depend. The two Eqs.(||) and (||) for pres- 
sure P must be equal for all saturations S. Combining 
the two gives a first order differential equation that dif- 
ferentiated once again with respect to S becomes Eq. ([!]). 
Note that in Eq. (Q) a = a' /A and b = b'/A. The second 
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FIG. 3: We show P(F nv ,) for a fixed value of the capillary 
number, namely Ca = 3.2 x 10~ 3 , which corresponds to Fig. 
^|(a). We study one sample assigning five different viscosity 
ratios M to the phases. We observe how it is possible to fit 
all the curves to the form in Eq. (|) for M € (1/30, 30). 




FIG. 4: The figure shows how four samples of fractional flow 
curves from our simulations can be fitted to the general solu- 
tion of Eq. (0) as it is given in Eq. (Q). The dimensionality, 
viscosity ratio and capillary number of these samples are as 
indicated in the legend. For clarity the curves are shifted so 
that the two curves in 2D belong the lower x-axis, and the 
two curves in 3D belong to the upper x-axis. The 3D curves 
are from the same simulation series as the 3D curves in Fig.nl. 



order version being more elegant is, however, autonomous 
and just as simple to solve using standard methods 
The general solution is 



F(S) 



1 (b-k) + (b + k) exp [k(S - Sp)} 
2a 1 + exp [k{S - S )\ ' 



where 



4afco, 



(4) 



(5) 



and fco and Sb are constants of integration. This result 
is generic and general. In Fig.^| we provide the fractional 
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FIG. 5: Data from the literature: Sharma and Yen(SY) p2fl, 
Peters and Khataniar(PK)[Q],Braun and Blackwell(BBj[g, 
p^ |. The water saturation is also the wetting saturation S w , 
but note that it is partly the different normalization of this 
entity, that was employed in the respective papers, which de- 
cide the placement of the curves. The SY-curves are typical 
samples of the set of functional forms that are employed when 
discussing how various physical parameters influence on frac- 
tional flow. Regarding PK and BB-curves, see Fig.^j. All the 
curves are fitted to the solution of Eq. (jll) as it is given in 
EqM. 




FIG. 6: The curves correspond to the respective curves of 
Fig.^. For PK the reported value of the viscosity ratio M — 
107.2 is used for the calculation from relative permeabilities 
to fractional flow and pressure. However, for BB we have 
chosen two different reasonable values for the viscosity ratio, 
to illuminate its effect. 



flow by simulations for four different sets of parameters. 
The sets are fitted by the form in Eq. (^) . We observe that 
the fits are excellent. Even though the two underlying 
observations are valid only in a certain central region 
(say 70-90%) of the curves, the solution is satisfactory in 
the entire two-phase region. 

Relative permeability is used both for steady-state flow 



and unsteady flow. The literature covers several meth- 
ods to produce relative permeability from experimental 
data||, 0|. It is a priori not entirely clear to which ex- 
tent a given study is comparable to our simulations. We 
have chosen a few samples from literature which are pro- 
vided in Figs.H and [| In literature, when focus is on 
how some physical or chemical properties influence on 
fractional flow, this is often done by employing a class 
of functional forms for the fractional flow. The curves 
indicated by SY are samples of that|2^]. It is reassur- 
ing that these curves can be fitted by Eq.^. The other 
three curves are calculated from the relative permeabil- 
ities fc mw and k rw and the viscosity ratio M using the 
formalism P — (fc rnw + Mfc rw ) _1 and F nw = Mk TW P, 
which arise when assuming that the total flux is held 
constant under the series and P is normalized pressure, 
for details see(|, 20 . 



The PK-curve is a steady-state curve which is com- 
parable to our simulations, even though the ensemble is 
constant total flux. The validity of the Eq.(p|) follows 
from Fig.|^, likewise Eq.(||) follows from Fig.^Here we 
use the actual viscosity ratio that was used in the ex- 
periments for our calculations. It is well established that 
relative permeability is a function of a large number of 
parameters including the viscosity ratioQ. However, to 
some extent one gets the impression that the curves are 
regarded as rock properties and are valid for at least a 
range of viscosity ratios. This might be very system de- 
pendent. The curves marked BB||, |2^] in Figs.|| and ^ 
are generated from relative permeability data in this way 
by choosing two values for viscosity ratio: M — 20 and 
M = 40. Within this range the results are robust. 

In general fractional flow curves are obtained for con- 
stant total flux. Another possible ensemble is constant 
pressure drop over a sample. In our simulations we have 
chosen a third ensemble, namely constant capillary num- 
ber. By this method we keep the ratio between viscous 
and capillary forces fixed. The precision and control of 
parameters we obtain in the curves exceed what is nor- 
mal in experiments. This is highly advantageous when 
looking for general relationships between variables. We 
believe that constant Ca is the appropriate ensemble for 
the study of steady flow properties. The comparisons 
with literature in this study indicate the correctness of 
the result, however, carefully designed laboratory exper- 
iments should be made using this ensemble to check our 
results. 

In conclusion it is our claim that this general result 
is valid for all steady-state two-phase flow system that 
can be modeled by our numerical model. That is to say 
immiscible flow where film flow can be neglected. This 
is more likely to be case at moderate to high capillary 
numbers. Studies of imbibition and drainage by small 
scale experiments show that in particular imbibition like 
steps may be dominated by film flow at low capillary 
numbers, but not at high capillary numbers p5[. 



5 



H.A.K. acknowledges support from VISTA, a collabo- 
ration between Statoil and The Norwegian Acad, of Sci- 
ence and Letters. This work has received support from 
NTNU through a grant of computing time on the good 
supercomputing facilities at NTNU. 



I — | * arendt@phys.ntnu.no 

1 Alex.Hansen@phys.ntnu.no 
[1] R. Lenormand, E. Touboul, and C. Zarcone, J. Fluid 

Mech. 189, 165 (1988). 
[2] D. H. Rothman, J Geophys. Res. 95, 8663 (1990). 
[3] T. A. Witten and L. M. Sander, Phys. Rev. Lett. 47, 
1400 (1981). 

[4] D. Wilkinson and J. F. Willemsen, J. Phys. A 16, 3365 
(1983). 

[5] F. A. L. Dullien, Porous Media: Fluid Transport and 
Pore Structure (Academic Press, San Diego, 1992). 

[6] D. Tiab and E. C. Donaldson, Petrophysics (Gulf Pub- 
lishing Company, Houston, 1996). 

[7] E. J. Peters and S. Khataniar, SPE Formation Evaluation 
pp. 469-474 (1987). 

[8] K. Langaas and P. Papatzacos, Transport in Porous Me- 
dia 45, 241 (2001). 

[9] M. S. Valavanides and A. C. Payatakes, Advances in Wa- 
ter Resources 24, 385 (2001). 
[10] M. Sahimi, Flow and Transport in Porous Media and 



Fractured Rock (VCH Verlagsgesellschaft mbH, Wein- 
heim, 1995). 

[11] M. J. Blunt, Current Opinion in Colloid & Interface Sci- 
ence 6, 197 (2001). 

[12] E. W. Washburn, Phys. Rev. 17, 273 (1921). 

[13] G. N. Constantinides and A. C. Payatakes, J. Colloid 
Interface Sci. 141, 486 (1991). 

[14] G. N. Constantinides and A. C. Payatakes, AIChE Jour- 
nal 42, 369 (1996). 

[15] J. Koplik and T. J. Lasseter, SPE J 25, 89 (1985). 

[16] E. Aker, K. J. Mal0y, and A. Hansen, Phys. Rev. E 58, 
2217 (1998). 

[17] E. Aker, K. J. Mal0y, A. Hansen, and G. G. Batrouni, 
Transport in Porous Media 32, 163 (1998). 

[18] H. A. Knudsen, E. Aker, and A. Hansen, Transport in 
Porous Media 47, 99 (2002). 

[19] S. Roux (unpublished). 

[20] H. A. Knudsen and A. Hansen, Phys. Rev. E (2002). 
[21] H. A. Knudsen and A. Hansen, in preparation for Phys. 
Rev. E. (2002). 

[22] M. M. Sharma and T. F. Yen, SPE Journal pp. 125-134 
(1983). 

[23] E. M. Braun and R. J. Blackwell, SPE-AIME 56th Con- 
ference, Texas (1981). 

[24] C. M. Bender and S. A. Orszag, Advanced mathemati- 
cal methods for scientists and engineers (McGraw-Hill, 
1978). 

[25] J.-D. Chen and J. Koplik, J. Colloid Interface Sci. 108, 
304 (1985). 



